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pq I We discuss two cases that can be connected to the dynamics of interacting populations: 

(I.) density waves for the case or neghgible random fluctuations of the populations 

^ ■ densities, and (II.) probability distributions connected to the model equations for of 

Y^- spatially averaged populations densities for the case of significant random fluctuations 

T^ ', of the independent quantity that can be associated with the population density. For 

c^ I the case (I.) we consider model equations containing polynomial nonlinearities. Such 

nonlinearities can arise as a consequence of interaction among the populations (for the 

case of large population densities) or as a result of a Taylor series expansion (for the case 

of small density of interacting populations) . In the both cases we can apply the modified 

t;;-!- ! method of simplest equation to obtain exact traveling-wave solutions connected to 

^ I migration of population members. Such solutions are obtained for systems consisting 

of 1 or 3 populations respectively. For the case (II.) we discuss model equations of 

.^ ■ the Fokker-Planck kind for the evolution of the statistical distributions of population 

O ■ densities. We derive several stationary distributions for the population density and 

^^ ' calculate the expected exit time connected to the extinction of the studied population. 

Keywords: interacting populations, nonlinear waves, method of simplest equation, 
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1. Introduction 

Since the famous paper of the properties of the Lorenz attractor half a century 
ago [l| the nonlinearities are intensively studied in different areas of science. Just 
several examples are (2| from the optics, [S] from biomechanics, [4] from solid state 
physics, |5|-[7| from fluid mechanics, |8|,|9| from the time series analysis, etc. Population 



dynamics is a classic area of application of nonlinear models [10|, lUj. In many cases 

the dynamics of interacting populations is studied on the basis of mathematical models 

consisting of equations that contain only time dependence of the population densities 

12|-(l7|. These models are very useful for understanding the complex djTiamics of 
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the interacting populations but they do not account for two important aspects of this 
dynamics: (I.) the possible influence of spatial characteristics of the environment; and 
(II.) the possible fluctuations of the population densities caused by different factors. 
Below we shall investigate two kinds of population dynamics models that account for 
each of these effects. First of all we shall discuss the dynamics of spatially distributed 



populations and this will be a continuation of our previous work [18|, [IQj . Then we shall 
show that by appropriate averaging the spatial model can be reduced to model in which 
the population densities depend only on the time. The models that contain spatially 
averaged quantities are valid for the cases of small and large values of the densities of 
the interacting populations. Because of this in the case of the model for the evolution 
of the spatially averaged population densities we shall discuss the influence of random 
fluctuations of the population densities for arbitrary values of the densities. The result 
of the action of these fluctuations is that instead of equations for the trajectories of 
the populations in the phase space of the population densities we shall write and solve 
equations for the probability density functions of the population densities. 

The organization of the paper is as follows. In Sect. 2 we obtain the model equa- 
tions for a system of interacting spatially distributed populations. In Sect. 3 after 
appropriate averaging the model system of nonlinear PD Es from Sect. 2 is reduced 
to a system of nonlinear ODEs for the temporal evolution of the densities of inter- 
acting populations. In Sect. 4 we obtain traveling- wave solutions of the system of 
model equations that can be associated with density population waves. In sect. 5 
we discuss the influence of fluctuations on the evolution of population densities of the 
model system of spatially averaged equations from Sect. 3. The inclusion of random 
fluctuations transforms the model system of deterministic nonlinear ODEs to a sys- 
tem of Langevin equations which is further transformed to a system of Fokker-Planck 
equations for the evolution of the probability density functions (p.d.f.s) of interacting 
populations. We discuss several examples for stationary p.d.f.s that are attractors for 
the time dependent p.d.f.s for the case of large times. Several expected exit times 
connected to extinction of populations are calculated too. In addition to the main 
text of the paper there are 3 appendices that are devoted to: (I) the modified method 
of simplest equation used to obtain the traveling-wave solutions; (II) description of a 
MAPLE program for obtaining the analytic form of the solution for the coupled kink 
waves for the case of 3 interacting populations; and (III) remarks on two observations 
connected to the diffusion Markov processes and the Fokker-Planck equation. 

2. Model equations 

2.1. Spatially distributed populations 

Let us consider an two-dimensional area S where A^ competing populations are 
present. The density of each population is pi{x,y,t) = ^^, where ANi is the number 
of the individuals of the i-th population that are present in the small area AS at the 
moment t. Now let a movement of population members through the borders of the 
area AS be possible and let ji{x, y, t) be the current of this movement. Then (j, ■ n)6l 
is the net number of members from i-th population, crossing a small border line SI 
with normal vector n. Let the density changes (other then these caused by the border 
crossings) be summarized by the function Ci{pi, p2, ■ ■ ■ , Pn,x, y, t). Then the change of 



the density of members of the i-th population in the studied area is described by the 
equation 

^ + divi = a. (2.1) 

Below we shall discuss the case where ji has the form of linear multicomponent diffusion. 
In this case 
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where Dik is the diffusion coefficient. Eq. (12.21) reflects the possibility that the mo- 
tion of the population members is caused not only by gradients of the density of the 
own population but also it could be caused by gradients of the densities of the other 
populations. 

We shall not specify the kind of the function Cj as we shall consider the general 
case of relatively small populations densities that allow us to write C, as Taylor series 
expansion around the zero values of all population densities as follows 
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We shall discuss the one- dimensional case and in addition we shall assume that the 
diffusion coefficients Di^ are constants. Then the substitution of Eqs. (12. 2p and (12. 3 p 
in (12. ip leads to the following system of nonlinear PDEs for the studied N interacting 
populations: 
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For the case of 1 population the system (12. 5p is reduced to the equation 
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Let us give two examples of systems of ki nd (12.511. The first example is connected to 
the classical Lotka- Volterra case extended in IJ] - 16| . In this case after assumption of 
dependence of the growth rates and cornpetition coefficients on the population density 



one arrives at the system of equations |14l|-|16 
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where i = 1,2, . . . ,n is a. number that indexes the n competing populations; Dij is the 
diffusion coefficient; A = -§^ + -§-^] rf and a°- are the parts of the corresponding growth 
rates and competition coefficients that do not depend on the population densities; Tj^ 
and aijk are parameters that regulate the intensity of the dependence of the population 
growth rates and competition coefficients on the population densities pi. It can be 
easily checked that the system described by Eq. (12. 7p is a particular case of the system 
of equations given by Eqs. fl2.ip - fl2.4l) . 

A system of kind similar to the system from Eq. (12.70 arises in the social dynamics 



in the spatial model of ideological struggle developed in 20| . For this case the model 
system of equations is 

rj n n n n 

-o^ - ^ DijApj = Tipi + J^ fijpj + J^ a^jPiPj + ^ bijkpipjpk + ■■■ (2.8) 
i=i j=o j=0 j,k=0 

where A = -^ + ^; pi,i = 1,2, ... ,n are the spatial densities of the populations of 
the followers of the corresponding ideology; r^ are the rates of change of corresponding 
populations of adepts by births and deaths; fij is the coefficient of non-contact conver- 
sion (the ideology of a person can be changed without contact between humans but by 
the mass media influence for an example); aij is the coefficient of binary contact con- 
version that describes the change of ideology by contacts between followers of different 
ideologies. 

2.2. Spatially averaged equations 

Below we shall apply spatial averaging to the system o f eq uations (12. 5p similar to 



the averaging used in the optimum theory of turbulence [21| - [23[. In the general 
two-dimensional case let a quantity q{x, y, t) be defined in an large two-dimensional 
plane area S with acreage \ S \. Then by definition the spatial average of q is 

qit) = -— / dx dy q{x,y,t). (2.9) 

q{x,y,t) can be separated in spatial averaged part q and the rest quantity Q{x,y,t): 

q{x,y,t)=q{t) + Q{x,y,t). (2.10) 

Let I S" I be large enough so that the plane average of any product of the rest quan- 
tities vanish: Qi = QiQj = QiQjQk = ■ ■ ■ = [J. In addition we shall assume that 
/ Jg dx dyV'^Q has finite and small value such that V^Q = jk J J^dx dyV^Q — )■ 0. 
The application of the averaging to Eq. (12. ip in presence of the assumptions given by 
Eqs. (12. 2p . (12. 3 p and (12. 4 p (note that in this case we have two spatial dimensions) leads 
to the system of ODEs as follows {i = 1,2, . . . , N): 

-, oo oo oo 
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^The products QiQj\ QiQjQk have finite values and by definition of tlie averages very large | 5' 
is presented in the denominator of the averaged quantity. 



We note that Eq. f l2.1ip follows directly from Eq. (I2.5p in the spatially homoge- 
neous case. The above discussion shows that Eq. fl2.1ip can arise also in the spatially 
inhomogeneous case. For the case of one population Eq. fl2.1ip becomes 
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We note that the equations of the kind of Eqs.f l2.ir]) and fl2.12p are often used as model 
equations in population dynamics not only for small values of population densities but 
also for large values of these densities, i.e., for large pj. One example is connected to 
Holling functional response functions in predator-prey systems [24|. For the case of one 
predator and one prey the functional response can be for an example type II Holling 

op 



functional response: /(p) 



p^ where p is the prey density and a and h are 
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parameters. The functional response can be also type III Holing functional response: 
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The functional response function for the case of many prey species can 



be more complicated. 

For small values of population densities the models even with complicated nonlinear 
functional responses can be reduced to models with polynomial nonlinearities. An 
example is the one predator - two prey model 25 
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where fi are the functional responses; m is the constant mortality rate of the predator; 
Cj are the conversion factors of captured prey species into predators and Qi are the 
growth functions of the corresponding prey type, p^ 2 ^^^ the spatial densities of the 
two types of prey and pg is the spatial density of the predator species. For small 
population densities we can apply Taylor series expansion for the functions /1.2 and 
gi^2 and thus obtain a system of equations from the studied in this paper class of 
equations 



dpi 
dt 

dp2 
dt 

dps 
dt 



/ . / . „ i„ 1 (Pl«3,ni,n2 -P3'^l,ni,«2) 



ni=0 712=0 



00 00 



ni\n2\ 



-ni—n2 



7 , 7 , 1„ 1 (P20^4,ni,n2 — P3'^'2,ni,n2) 



ni=0 712=0 
P3 



ni\n2\ 



00 00 — m — 712 
Pi P2 



-m+J2 Yl 



ni=0 712=0 



ni\n2\ 



{ClOil,ni,n2 + C2«2,ni,7i2) 



(2.14) 



where 
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3. Traveling waves 



3.1. Case of one population 

Let us discuss the simplest case of one population described by Eq. (12. 6p . First we 
introduce the travehng-wave coordinate ^ = x — vt where v is the velocity of the wave. 
In addition we shall assume that the polynomial nonlinearity in Eq. (12. 6p is up to order 
L. We rescale the coefficients in Eq. (l2.6p as follows: 

D'^ = -D/v] al^ = anjv. 

Then Eq. (l2.6p becomes: 
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Below we shall obtain exact solution of Eq. (13. 2p by application of the modified method 
of simplest equation for obtaining exact solutions of nonlinear PDEs. For more details 
on the modified method of simplest equation see Appendix A. 

Proposition 1. The balance equation for Eq. li3.S\) for the case when Riccati equation 
is used as simplest equation is P{L — 1) = 2 where P is the largest power in the 
polynomial for p[^) constructed on the basis of the solutions $(0 of the Riccati equation. 

Proof. We apply the methodology from Appendix A to Eq. (13. 2p . We constrict a solu- 
tion as finite series 
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where $(^) is a solution of the Riccati equation 
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arise from the different terms of Eq. (l3.2p (these powers are P + 2 from the term ^ 



The substitution of Eq. (l3.3p in Eq. (l3.2p and the balance of the largest powers of $ that 

d^J) 
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and PL from the term a\p^) lead to the balance equation 

P(L-1) = 2. (3.6) 
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Thus we have the possibihties: P = L 
possibihties. 



2 or P = 1: L = 3. Below we discuss these 



3.2. Case P = L = 2 

In this case we shall obtain exact traveling-wave solution of the equation 
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Proposition 2. The traveling wave solution of the kind Eq. li3.3\) of Eq. ( [J. 7[ j obtained 
on the basis of the modified method of simplest equation when the equation of Riccati 
is used as simplest equation is 
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for the case aj 7^ and 
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/or i/ie case aj = 0. .^0 is a constant of integration. 

Proof. As P = L = 2 then from Eq. (13. 3p the solution will be of the kind 
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The substitution of Eq. (l3.10p in Eq. (l3.7p leads to a system of relationships for the 
parameters of the solution (the system is of kind (lA.4p ). 
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Now we have 2 possibilities: aj 7^ and aj = (which is more close to the classical 
population dynamics models that usually do not possess terms independent on the 
population density). 

Case aj 7^ 

For this case the solution of the system (13. lip is as follows: 
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and then we obtain the solution fl3.8p of Eq. (13. 7p 
Case Oq = 

For this case the solution of the system (13. lip is as follows: 

6 3662 _^ gQ^^t _^ 25af 



(3.12) 



D^ 



Oi 



25a\ ' 
(3662 



flo 



100a\al 



25a[ )(66 + 5al 



Q00ca\a 



t^.t 
2 



02 



(3662 _ 25c^t ^ 
14400cala^ 



3662 _ 25c^t 

a = . 

144c 

and then we obtain the solution (13. 9p of Eq. (l3.7p . 
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The obtained solutions describe kink waves that can be considered as traveling 
waves of change of the value of the population density of the studied population. Such 
waves can describe the front of migration of a population. The appropriate values of 
the boundary conditions ensure that p{C,) is non-negative elsewhere. We note that the 
parameters of the solved Eq. (l3.7p are D^^ and aQ,a\,a2- The first relationship from 
Eq. ( 13.12P connects these 4 parameters. Then the solution ( 13. Sp does not hold for any 
values of parameters of Eq. (l3.7p but only for these combinations that satisfy the above 
mentioned relationship. 

Let us discuss in several words the problem with the boundary conditions of the 
obtained solutions, theoretically for the general case of solution (13. 8p there are 10 
parameters and 5 relationships (13. lip among them. Thus there are 5 free parameters 
and several possibilities for boundary conditions are 
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Figure 1: Several solutions that satisfy boundary conditions p{+oo) = Ai; p(— oo) = A2. Solid 
line: Ai = 1.5; A2 = 0.5. Dashed line: Ai = 1; A2 = 0. Dot-dashed line: Ai = 1.5; A2 = 0. 
Dot-double-dashed line: Ai = 0; A2 = 2. 



Let us impose the boundary conditions p(-foo) = Ai; p(— 00) = A2 on the solution 
given by Eq.f l3.8p . The result is that there are two additional relationships that must 
be satisfied by the parameters of the solution. The relationships are as follows 
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The second relationship can be written as follows: A2 = —{Ai + ayal) which means 
that if the parameters a| 2 of the solved PDE are fixed then when we choose the 
boundary condition Ai then the boundary condition A2 can not be arbitrary. For an 
example ii Ai = then A2 = —al/al- Several examples of the obtained nonlinear 
waves satisfying two boundary conditions are shown in Fig.l. 



3.2.1. CaseP = 1; L = 3 

In this case we shall obtain exact traveling-wave solutions of the equation 
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Proposition 3. The traveling wave solution of the kind Eq. \3.'J\) of Eq. \3.1b]) obtained 
on the basis of the modified method of simplest equation when the equation of Riccati 
is used as simplest equation is 
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Proof. As P = 1;L = 3 then from Eq. (l3.3p the solution will be of the kind Eq. (l3.17|) 
The substitution of Eq. (l3.17p in Eq. (l3.16p leads to the following system of relationships 
for the parameters of the solution 
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We shall consider two possibilities. First we shall solve the system (I3.2ip with respect 
to q:qi2 3. This means that the coefficients ao,i as well as the coefficients a, b, c of the 
Riccati equation will remain free and we can impose many boundary conditions on 
the solution. Thus we shall investigate a small subclass of equation of kind Eq. (l3.16p 
where we have solution with many free parameters. Then we shall consider the case in 
which we shall keep as many as possible of the parameters of Eq. (l3.16l) free. The price 
for this will be that we shall have to impose several relationships on the coefficients of 
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the solutions 

Case 1: Solution of 113. 2 1\) with respect to al ^23 

This solution is as follows 
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Thus we arrive at solution f l3.17p of Eq. f l3.18p . 

Case 2: Solution of Ii3.21\) with respect to a, oq, ai and aj 

The solution is as follows 
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Thus we arrive at solution (13.191) of Eq. (I3.20p . D 



We can now impose boundary conditions of the solution (13.171) of Eq. (I3.18P . For 
an example the boundary conditions can be 



p(+oo 

dp 



Ai] p(-oo) = 
^3, (5i > 0); 



B^ 



^2; 

dp 
d^ 



A,, {B2 > 0) 



-B2 



The boundary conditions will fix additional parameters of the solution. We let the 
corresponding algebraic manipulations to the interested reader. 
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3. 3. Coupled waves in a system of 3 populations 

Let us discuss a system of 3 competing populations modeled by the system of 
Lotka-Volterra kind f l2.7p for the case of constant coefficients of change of population 
members and constant interaction coefficients, i.e., for the case Tj^ = 0; aijk = 0. The 
system of equations becomes 



^1 n ^1 



dt 
dt 
dt 



D21 
D31 



dx 

dx 
^1 



ox ox 

i^22-^ i^23-^r- 

OX OX 



r\px - r\a\^p\ - r\ot\^p^p^ - r\a\^p^p.^ 



AP2 - rl(^\\P\p2 - A(A2P\ - r2^23p2p3 



r^ ^P2 p. ^PS _ 0- ^0^0 -^ -j: „0^0 -j: -j: „0^0 -j:2 

dx "aJT ~ ~dx ~ ^^ ~ 3"3lPlP3 - ^3"32P2P3 " ^3"33P3 



(3.24) 



For this system we shall demonstrate existence of a simple coupled kink wave solution 
(more complicated solutions are possible too). 



Proposition 4. The system i3.24\) possesses coupled kind wave solution of the kind 



Pi(0 

P2(0 
P3(0 



«iCo I ai + 61 + 4aiCo 9(ai + 61 + 4aiCo) ^ , 
+ cti < — : — : ^- \ : — — ^^ — tanh 



ai + hi 
bicp 



+ &1 



4ai(ai + 61) 

ai + bi + 4aiCo 

4ai(ai + hi) 



4ai6(ai + hi) 
(ai + 61 + 4aiCo) 
4ai6(ai + 61) 



tanh 



, fli + &i + 4aiCo 6^(01 + 61 + 4aiCo, , , 
-{ai + 61) + ci <; — ^ — 7 — p^ 1 ; — — r^ — tanh 



2 

' ^(e + ^o) 
2 

^(e + eo 



where $, — x + 



4ai(ai + 61' 



4aiCo + (ai + 6i)(l - 6 + fo^n; 
b{ai + 61) 



4ai6(ai + 61 



(3.25) 



■t. fli, 6, hi, Co, -Dii are free parameters. 



We remember that above 6"^ = b^ — Aac where a, 6, c are the parameters of the 
Riccati equation. 



Proof. We apply the modified method of simplest equation to the system (I3.24p . First 
of all we introduce the traveling wave coordinate ^ = x — vt where v is the wave velocity. 
Then we search for the solution in the form 

p Q R 

MO = J2^^H0'■, MO = J^^MO'; MO = J^^kHO' (3.26) 



i=0 



fc=0 



d^ 

where — — = a$ + &$ + c. The simplest possible balance equation is P = Q = R = 1. 

The substitution of all above in the system (I3.24p leads to the following system of 9 
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5.) 



nonlinear algebraic equations 

1-) --DisCia - {v + Du)aia - noauoaibi + rioanofli - Duha - rioaisoaiCi = 
2.) -rioauoaobi - rio«i2oai&o + 2rioaiioaoai - rioCHsoaoCi - rioai3o«iCo - no«i - 

Discib - Dubib -{v + Dn)aib = 
3.) -{v + Du)aic - riQai2oaobo - -D13C1C - D^bic - rioaisofloCo - ^loao + 

rwauoal = 
4.) -D23Cia - D2iaia - r2oa220&? + r2o«2ioai&i - {v + D22)bia - r2oa23obiCi = 

-2r2oa22obobi + r2oa2ioao&i + ^2oa2] 

-D23C1& - (t; + D22)bib - D2iaib = 
6.) -D2iaic - r2oa22obl - D23C1C - {v + D22)bic - r2oa230&oCo - r2o&o + r2oa2ioao&o = 
7.) -(t; + D33)cia - D-^aia - r^oa^iobiCi + rsoasiocuci - /^32&ia - r^oa^^ocl = 
8.) -r3oa32o&oCi - r3oa32o&iCo + r3oa3ioaoCi + r3oa3ioaiCo - 2r3oa33oCoCi - tsqCi - 

{v + D33)cib - -D32&1& - Dsiaib = 
9.) -Dsiaic - r3oa32oboCo - {v + £'33)cic - £'32&ic - r3oa33oCo - tsqCo + r3oa3ioaoCo = 

(3.27) 

The general solution of this system is very long. In order to obtain the solution from 
the text of Proposition 4 we fix some of the parameters in the above system as follows 



-D23Cia - D2iaia - r2oa220&i + r2o«2ioai&i - (v + D22)bia - r2oa23obiCi = 1 
-2r2oa22obobi + r2oa2ioao&i + ^2oa2ioai&o - '^2o"230&oCi - r2oa230&iCo - ^20^1 

DnoCA b — (iJ -\- Dnn)bi b — Dm nib ^ Q 



no = i;'^2o = i;^3o = i;«iio = i;«22o = i;a33o = i;«i2o = i;«i3o = i;«2io = 1; 

^230 = 1; ^310 = 1; ^320 = 1; D21 = D12; -D31 = -D13; -D32 = -D23; -D22 = Du] D33 = Du] 

D,2 = l;Dr3 = l;D23 = l. 

(3.28) 

In this simple case the solution of the system (I3.27P is (for details of solution of the 
system of equations by means of a Maple program see Appendix B) 

ci = -(ai + 61) 

4aiCo + bDiiai — aib + ai — bib + bDubi + bi 



V 



Oo 



b{ai + bi] 

Co(2aiCo + ai + bi)b 

(fli + 61) (ai + 61 + 4aiCo) 

aib{ai + bi) 

(ai + 61 + 4aiCo) 
aico 



ai + b 



bo = -^ (3.29) 

ai + 61 

Substituting the coefficients in the functions p7, p2, P3 we obtain the coupled kink wave 
solution from the formulation of the Proposition 4. D 
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4. Statistical distributions and exit time 



Eqs. fl2.lip and f l2.12p are typical equations for description of the evolution of dy- 
namical systems. The general case of such equations is 



dt 



XAxi, X2, 



^Xn 



1,2, 



A^. 



(4.1) 



where Xj(xi, X2, ■ ■ ■ , xn) is some (in the general case nonlinear) function. For such kind 
of systems there exists a theory that allows us to characterize some system properties 
in the case when the system is under the action of random perturbations. Pontryagin, 
Andronov and Vitt |4]J] developed such theory for random impulses that occur after 
every interval of time r and each impulse causes the phase point of the dynamical 
system described by Eqs. (14. ip to jump through a distance a along a random direction. 
Let us first consider the case of single population and one spatial dimension. For the 
case when a tends to together with r in such a way that the ratio a^/r tends to finite 
limit b it is possible to obtain an equation for the probability density function p{x, t) 
as follows (for more discussion see Appendix C): 



dp d 

dt dx 



X{x)p 



b d'^p 
2dx^' 



(4.2) 



For the general case of A^ populations the equation for the probability density function 
becomes 



N 



dp ^^ d 



dt 



i=\ 



dxi 



Xi{xi,X2,...,XN)p 



N N 



i=l j=l 



d'^p 

dxidxj ' 



(4.3) 



where bij are again coefficients that characterize the random impulses. 

Another kind of problem that can be solved by this approach is to calculate the 
mathematical expectation of the exit time. Let us again first discuss the case of one 
population and one spatial dimension. We have a phase point that is inside the interval 
[ei, €2] (ei < €2) and the system is under the infiuence of the same random perturbations 
as described above. The exit time is the time for which the phase point that was inside 
the above interval at t = will leave this interval through ei or through €2- If we 
denote as F{x) the mathematical expectation for the exit time then F{x) is a solution 



of the equation |41l|.|42 



b d^F dF ^ 

2 dx'^ dx 



0, 



(4.4) 



with boundary conditions F{ei) = -^(62) = 0. For the case of many populations the 
zero boundary conditions are on the entire border of the multidimensional phase space 
area that has to be exited and the equation for the probability density function of the 
exit time is 



N N 






d'^F 

' dxjdx. 



N 



+ ^X(a;i,X2, 



i=l 



,dF 



(4.5) 



Let us now apply this theory to Eq. ( 12.121) . We shall be interested in the stationary 
distributions p{j>) connected to Eq. (j2.12p .i.e.. we shall be interested in the case when 
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after a long time the probability density function p becomes stationary and depends 
only on the spatial coordinate p. This stationary case is important because of Obser- 
vation 2 from Appendix C which states that each time dependent solution p{x, t) of 
the Fokker-Planck equation (1C.7|) converges at t — )■ oo to the stationary distribution 
Po{x) from Eq. (IC.12| . In our case x = p and X{x) = ^^ =o'^niP"^ We shall discuss 
two cases: (i) the case of single population where p G [0, oo); and (ii) another case not 
connected directly to the population dynamics where p G (—00, 00). 
Case 1: p E [0, 00) 

This case is connected directly to the population dynamics as the population density 
can't be negative, in this case we have 



Proposition 5. Let us discuss a system described by the equation 



(4.6) 



ni=0 



which is under the action of random impulses that occur after every interval of time 
T and each impulse causes the phase point of the dynamical system described by Eqs. 
( [^.i^ to jump through a distance a along a random direction. Let when a tend to 
together with r in such a way that the ratio a?' jr tends to finite limit b. Let in addition 
the following requirements be fulfilled 



Pi0) = 0; pG [0,00) 
Then the stationary p.d.f. p(p) is 



p{p) = Cexp 



E 



ni+l 



a 



"1 



ni=0 



ni + 1 



L 



2 >:^ p"i + l 

ni + 1 



ni=0 



(4.7) 



(4. 



dp exp 

where the constant of integration C is determined by the normalization condition 

dp p{p) = 1. (4.9) 



Proof. The equation of the stationary distribution p{p) is particular case of Eq. (14. 2p . 
One integration of the equation for the stationary distribution leads to 



b dp 
2lp 



Pip) Yl ""1^"' 



ni=0 



Ci=0 



(4.10) 



Eq. (I4.10p can be easily integrated and the result is 

, L 



pip) = exp 



XE""! 



■Tflll + l 



■pn. 



ni=0 



rii + 1 



c 



2Ci 



dp exp 



h Z^""i 



Tr"i+1 



■pni 



ni=0 



Tli + l 



(4.11) 



where C is a constant of integration. The condition p{0) = leads to Ci = bC/2 and 
distribution described by Eq. (14. lip is reduced to the distribution from Eq. (14. 8p . D 
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Case 2: p G (— C)0, oo) 
Proposition 6. Let us discuss a system described by the equation 

L 



dp 



i^E"...?" 



(4.12) 



ni=0 



which is under the action of random impulses that occur after every interval of time 
T and each impulse causes the phase point of the dynamical system described by Eqs. 
^■1^ ) to jump through a distance a along a random direction. Let when a tend to 
together with t in such a way that the ratio a^/r tends to finite limit b. 



* Subcase 1: 

Let in addition the following requirements be fulfilled 

p{0) = A; p E (— oo, oo) 

Then the stationary p.d.f. p(p) is 



(4.13) 



p{p) = exp 



h Z^ ""1 



7T"1 + 1 



■pni 



ni=0 



ni + 1 



' h Z^ ""1 



TT^l+l 



pni 



ni=0 



rii + l 



C -{C-A) j d-p exp 

where the constant of integration C is determined by the normalization condition 

dpp{-p) = l. (4.15) 



4.14) 



irk Subcase 2: 

Let in addition the following requirements be fulfilled 

1 dp , 2ao _ / X 

^ |p=o= ^-; pe (-00,00) 



p(0) dp '"-" b 

Then the stationary p.d.f. p(jj) is 



(4.16) 



p{p) = C exp 



L 



2 -^ p"i+^ 



h ^ "' ni + l 

ni=0 



(4.17) 



where the constant of integration C is determined by the normalization condition 

dpp{-p) = l. (4.18) 



Proof. Subcase 1 

The equation of the stationary distribution p(p) is particular case of Eq. (14. 2p . One 
integration of the equation for the stationary distribution leads to Eq. (I4.10p . Eq. 
(I4.iup can be easily integrated and the result is Eq. (14. lip , where C is a constant 
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P0.2- 




P0.2- 




Figure 2: Several profiles oi p{p) from Eg. (14.171) (in the figures p is denoted as p). Figure (a); b — 2; 
ao = 0; ai = 1; a2 = —0.4; a^ = 0.2; a^ = 0; as = —0.5. Figure (b): ao = 0; ai = 1; a2 — 0; 
as = 0.2; a4 = 0; as = —0.5; 6 = 2. Here the fixed points are 3. The two maxima of the p.d.f. 
distribution are centered around the two stable fixed points and the minimum is centered on the 
unstable fixed point p = 0. Figure (c): ao = 0.003564; ai = -0.006084; aa = 0.3975;a3 = -1.234; 
ai = 1.81; as = -1; 6 = 2. Figure (d): ao = 1.2936; ai = -7.2436; a2 = 14.58; ag = -13.63; ai = 6; 
as — —1; 6 = 2. The probability distribution has 3 maxima (at p = 0.4, 1.1, 2.1) and two minima. 



of integration. The condition p(0) = A leads to Ci = b{C — A)/2 and distribution 

described by Eq. (14. lip is reduced to the distribution from Eq. (I4.14p . 

Subcase 2 

The integration of Eq. (l4.2p leads to Eq. (l4.10]) . where Ci is a constant of integration. 

This constant is equal to because of the condition (I4.16p . With Ci = we can 

continue the integration of Eq. (l4.10p and the result is (I4.17P where the constant of 

integration C is determined by the normalization condition (I4.18p . D 



We note that Eq. (14.161) in combination with Eq. (I4.17P mean that ^„ ^g a„^ = 0. 
In addition f{p) must tend to when p — )• ±oo. The dominant term at large values 
of p is aip^. Then a^ must be negative (to ensure / — )■ at large positive values of p 
and L must be odd (to ensure / — )• at large negative values of p). 

Several examples for statistical distributions /(p) are shown in Fig. 2. Let us note 
here the tunnel phenomenon which arises because of the presence of fluctuations. The 
existence of the distributions p{p, t) and po{p) means that in the course of the time each 
value of the density p can be reached. Then if for an example at the initial moment the 
system trajectory in the phase space is close to a fixed point of the ODE without added 
fiuctuations then in the course of the time the phase point can leave this area and to 
travel to phase space area that is close to another fixed point. This tunnel phenomena 
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4 6 8 10 

P 




Figure 3: Exit time expectations for q = (which means extinction of the population) calculated 
on the basis of Ea. (l4.19p . p on the horizontal axis is equal to p from Eq. (|4.19p . For all curves 
6=2. Figure (a): influence of the value of a^ on the exit time, a^ — 0.01; ai — 0.2; 0:2 = 0.1. 
Solid line: a^ = —0.05; dashed hue: a^ = —0.04; dot-dashed line: a^ = —0.03; dot-double dashed 
line: a^ = —0.02. Figure (b): influence of the value of ai on th exit time, ao — 0.01;a2 = 0.1; 
a — 3 = —0.02. Solid line: ai — —0.3; dashed line: ai — —0.4; dot-dashed line: ai — —0.5; dot-double 
dashed line: ai — —0.6. 

is closely connected to the role played by fluctuations in the case when a bifurcation 
happens in the studied system. 

Let us now calculate the exit time expectation on the basis of Eq. fl4.4p . 



Proposition 7. let us discuss the system described by Eq. ^.l2^ . The distribution 
Fqij)) from the initial position p to the position q KJ) is 



F,ip) 



d^ exp 



?ii=0 



e 



ni+l 



ni + l 



dr] exp - ^ 



a„ 



ni=0 



-.rii+l 



rii + l 



(4.19) 



when 

n F{p = q) = , 

{**) -F(p, q) increases in the slowest possible manner as p —)■ 00 . 

Proof. We calculate the distribution for exit time from the initial position p to a po- 



sition q < p. In the discussed case again X{x) 
Eq. fl4.4p leads to the equation 



l^ni=0 



ttniP"'. 



One integration of 



dF 

-^ = exp( 

dp 



-i'{-p)){ci + 



d^-exp{^{0) 



^^^^-bj. 



d^XiO (4.20) 



The relationship (**) requires Ci = (as the corresponding term in Eq. fl4.20p vanishes 
and the growth of ^ is as slow as possible) , and the integration of Eq. (I4.20p leads to 
the result KWi D 

Fig. 3 shows the dependence of the exit time expectation on the population density 
and coefficients of the model equation for the case L = 3. The negative values of 
q;i,3 make extinction be expected sooner whereas the positive values of the other two 
parameters can delay the extinction. The theory can be easily applied for the case of 
system of many interacting populations but even in the simplest one-dimensional case 
the integral from Eq. (l4.19p must be calculated numerically. 
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5. Conclusion 

In this paper we liave discussed two aspects of population dynamics. First we have 
presented a model of the space-time dynamics of the system of interacting population 
in 2 spatial dimensions. For the most simple case of one spatial dimension and for 
one population we have obtained exact traveling-wave solution of the model nonlinear 
PDE by means of the recently developed modified method of simplest equation for 
obtaining exact and approximate solutions of nonlinear PDEs. The obtained exact 
solution describes the propagation of changes of the population density in the space. 
The generalization of this theory to the case of many populations is straightforward 
and describes the spreading of coupled waves of changes of densities of the studied pop- 
ulations. The case of three populations is discussed in the paper. The second discussed 
aspect of the population dynamics was connected to the influence of the random fluctu- 
ations on the population densities. The presence of fluctuations leads to description in 
terms of probability density functions for the population densities. The discussed gen- 
eral theory is illustrated again for the simplest possible case of one population in two 
aspects: calculation of probability density functions and calculation of the expected 
extinction time. As expected from the theory of diffusion Markov processes the min- 
ima and maxima of the obtained probability density functions are exactly at the fixed 
points of the corresponding non-perturbed model system of differential equations. The 
expected extinction time strongly depends on the coefficients of the model equations. 
Finally several results are obtained that are not connected to the theory of interacting 
populations but may be interesting for other cases modeled by nonlinear PDEs with 
polynomial nonlinearities. 

Acknowledgment 

This research was partially supported by the Fund of Scientific Researches of Re- 
public of Bulgaria under contract DO 02-338/22.12.2008 in the scope of which the 
averaging applied in section II B has been developed and used. 

Appendix A. Modified method of simplest equation 

There are many approaches for obtaining exact analytic solutions of nonlinear par- 



tial differential equations J26|- J32| . In this paper we use the modified method of simplest 



equation. The schema of application of the modified method of simplest equation is 



shown in Fig.4. The method of simplest equation [33|-[35j is based on the fact that 
after application of appropriate ansatz large class of NPDEs can be reduced to ODEs 
of the kind 

and for some equations of the kind (1A.1I) particular solutions can be obtained which 
are finite series 



F{o = J2^^mo]\ (A.2) 



j=0 



19 



Solved nonlinear PDE 


i 


Reduction ansatz to a 
nonlinear ODE 


'' 



Construction of a 
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Figure A. 4: The modified method of simplest equation 



constructed by solution $(.^) of more simple equation referred to as simplest equation. 
The simplest equation can be the equation of Bernoulli, equation of Riccati, etc. The 
substitution of Eq. (IA.2p in Eq. flA.ip leads to the polynomial equation 



•P = (To + 0-1$ + 0-2"^^ H \- O^r^' = 0, 



(A.3) 



where the coefficients a,- 



0, 1, . . . , r depend on the parameters of the equation and 



on the parameters of the solutions. Equating all these coefficients to 0, i.e., by setting 



(Ji 



0,^ = 1,2, 



r. 



(A.4) 



one obtains a system of nonlinear algebraic equations. Each solution of this system 
leads to a solution of kind (1A.2|) of Eq. ( lA.ip . 

In order to ensure non-trivial solution by the above method we have to ensure that 
ar contains at least two terms. To do this we have to balance the highest powers of 
that are obtained from the different terms of the solved equation of kind (lA.ip . As a 
result of this we obtain an additional equation between some of the parameters of the 



equation and the solution. This equation is called balance equation [36[-|40 . 
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Appendix B. Maple program for solving the coupled waves case 

Here we present a Maple program for obtaining exact solution of a system of equa- 
tions for description of coupled waves in a system of 3 populations discussed in Sec... 
The program has the following parts. 



Part 1: The equations 

l.> eql := —{v+Dll)*dif f(rhol(xi), xi)—D12*dif f{rho2(xi), xi)—Dl3*dif f(rho3{xi) 

rlO * rhol{xi) + rlO * alphallO * rhol{xiY — rlO * alphal20 * rhol{xi) * rho2{xi) — 

rlO * alphalSO * rhol(xi) *rho3{xi); 2.> eq2 := —D21*diff(rhoi(xi),xi) — (v + D22) * 

di f f {rho2{xi) , xi) — D23 * diff{rho3{xi),xi) — 

r20 * rho2{xi) + r20 * alpha210 * rhol{xi) * rho2{xi) — r20 * alpha220 * rho2{xi)'^ — 

r20 * alpha230 * rho2{xi) * rho3{xi); 

3.> eq3 := —D31*dif f{rhol{xi), xi)—D32*dif f{rho2{xi), xi) — {v+D33)*dif f{rho3{xi) 

r30*rho3{xi)+r30*alpha310*rhol{xi)*rho3{xi)—r30*alpha320*rho2{xi)*rho3{xi) — 

r30 * alpha330 * rho3{xiY; 



XI] 



XI) — 



Part 2: Substitution of the relationships for pi,2,3 in the equations 

4.> rhol{xi) := aO + al * Phi{xi); 
5.> rho2{xi) := 60 + 61 * Phi{xi); 
6.> rho3{xi) := cO + cl * Phi{xi); 
7.> eql] eq2; eq3; 



Part 3: Substitution of the relationship for ^ in the equations 



8.> eqla := subs{diff{Phi{xi),xi) ■ 
9.> eq2a := subs{diff{Phi{xi),xi) ■ 
10. > eq3a := subs{diff{Phi{xi),xi) 



a * Phi{xiY + h* Phi{xi) + c, eql); 
a * Phi{xif' + 6 * Phi{xi) + c, eg2); 

\2 



a * 



Phi{xi) + b* Phi{xi) + c, eq3); 



15.> e2 
16.>e3 



18.>e5 
19.>e6 



21.>e8 
22. > e9 



Part 4: Extracting the system of nonlinear algebraic equations 

11. > eqlb := collect{eqla, Phi{xi)); 
12. > eq2b := collect{eq2a, Phi{xi)); 
13. > eq3b := collect{eq3a, Phi{xi)); 
14.> el := coef f{eqlb, Phi{xi)'^); 

coef f{eqlb, Phi{xi)); 

— (f + Dll) * al* c — rlO * alphal20 * aO * 60 — D13 * cl* c — D12 *bl* c- 
rlO * alphal30 * aO * cO — rlO * aO + rlO * alphallO * aO^; 
17.> e4 := coeff{eq2b,Phi{xify, 

coef f{eq2b, Phi{xi))] 

-D21 *al*c-r20* alpha220 * 60^ - D23 *cl*c-{v + D22) *bl*c- 
r20 * alpha230 *bO*cO- r20 * 60 + r20 * alpha210 * aO * 60; 
20.> e7 := coeff{eq3b,Phi{xify, 

coef f{eq3b, Phi{xi)); 

—D31 * al* c — r30 * alpha320 * 60 * cO — (f + 1^33) * cl* c — D32 *bl* c- 



r30 * alpha330 * cO — r30 * cO + r30 * alpha310 * aO * cO; 
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Part 5: Simplifying assumptions 

23. > rlO := 1; r20 := 1; r30 := 1; alphallO := 1; alpha220 := 1; alpha330 := 1; 
alphal20 := 1; alphalSO := 1; alpha210 := 1; alpha230 := 1; alphaSlO := 1; 
alpha320 := 1; D21 := L)12; D31 := D13; D32 := D23; D22 := Dll; D33 := Dll; 
D12 := l;D13 := I; D23 := 1; 

24.> el; e2; e3; e4; e5; e6; e7; e8; e9; 

Part 6: Solution of the system of nonlinear algebraic equations 

25. > so/1 := solve{el,Ci); 

26.> cl := so/1; 

27.> so/2 := solve{e2,v); 

28.> V := so/2; 

29. > so/3 := so/fe(e3,c); 

30. > c := so/3; 

31. > so/4 := so/fe(e4,a); 

32.> a := so/4[l]; 

33. > so/5 := so/fe(e5, Oo); 

34. > flo '■= so/5; 

35. > so/6 := solve{e8,bo); 

36.> 60 := so/6; 

By this program one obtains relationships for the parameters aQ,bQ,Ci,a,c,v. 

Appendix C. Fluctuations, diffusion Markov process and Fokker-Planck 
equation 

In the paper we discuss equations of the kind 

- = X[x(t)]+r^(t) (Cl) 

where the process ri{t) models the small fluctuations. Let rjit) = (jC,{t) where a > 
is the intensity factor and the covariance function of the process ^(t) be a (5-function 
E[^it)C{t)] = S{t - t'). If in addition the expected value of Hf) is 0: E[i{t)] = the 
process ^(t) is called white noise and the equation 

— =X[xit)]+aaty, a;(0)=a;o (C.2) 

is called Langevin equation. 

^{t) can be written as time derivative of a Wiener process Wt (for this process 
Wo = 0; the function Wt(t) is almost surely continuous; and Wt has independent 
increments Wt — Wg (0 < s < t) which are normally distributed with expected value 
and variance equal to t — s): 

m = ^ ^Wt = J^ds^is) (C.3) 

Then the Langevin equation can be written in the form 

dxt = X{xt)dt + adWt] xq : random (C.4) 
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After one integration of the Langevin equation one obtains 

xtioj) = xo{u) + / dsX[xs{uj)] + aWtiu) (C.5) 



as the white noise is (5-correlated the solution of Eq. (IC.5l) is a homogeneous Markov pro- 



cess. The infinitesimal generator A of the solution process Xt |42|, [43| is a differential 
operator of second order 

(Ag)ix) = fix)g'ix) + ^g"ix) (C.6) 

The form of the infinitesimal operator A is important as it is determines the form of 
the Fokker-Planck equation for the one dimensional distribution p{x, t) connected to 
the solution process (IC.Sp . This Fokker-Planck equation is 

*fc^ ^ - A |X(.M.. 01 + ^%^; p(., 0) ^ p„(.) (C.7) 

We note that if we set a^ = 6 we obtain Eq. (14.21) from the main text of the paper. By 
setting 

J,=PX-!^| (C.8) 

we can write Eq. (1C.7l) in the form 



Eq. (lC.9P is a balance equation for the probability. There is no production of probability 



and the balance equation shows that the transport of probability along x happens by 

a"^ dp 
'flow', determined by the 'velocity' field X and 'diffusion' determined by 7—. This 

2 ox 
'diffusion' tries to decreases the differences in the probability distribution. 

Because of the above-mentioned the considered class of Markov processes are called 
diffusion processes. X is the drift of the diffusion process and a^ is called diffusion 
coefficient of the diffusion process. 

We finish this short discussion by two observations that are important for the text 



in the body of the paper [42 

Observation 1: 

In the stationary case p{x,t) — )■ p{x,0) = Po{x) the stochastic system described by 
Eq. ( [C.^I J spends much time around the stable fixed points of the corresponding deter- 



ministic system 

Observation 1 is a consequence of the fact that for the stationary case the Fokker- 
Planck equation (1C.7|) becomes 
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After one integration one obtains 



(T^ dpo{x) 



4 = -X{x)po{x) . 2 ^^ 



const 



(C.ll) 



The integration constant is equal to because of the boundary condition J^ = |±oo 
and for po{x) one obtains 



exp 



Po(x) 



a^ 



,V{x) 



J^ dx exp 



2 

—^V{x) 



o" 



(C.12) 



where 



V{x) 



dz Xiz) 



The local extrema of po(^) are given by 
dpo , 



dx 



\xo~ 



wrm dv 

dx 



\xo~ 



y f{xo) 







(C.13) 



(C.14) 



Thus the extrema oi po{x) coincide to the fixed points of the corresponding determin- 
istic differential equation. The maxima of Po{x) are located at the stable fixed points 
and the minima of ^0(2^) are located at the unstable fixed points. 

The second observation is connected to the importance of the stationary distribution 

Po{x): 

Observation 2: 



Each time dependent solution p{x,t) of the Fokker-Planck equation ( [C. 7| j converges at 
i — 7- 00 to the stationary distribution po{x) from Eq. /lC.l^) (if po{x) exists). 

In order to show that Observation 2 is true we shall use the technique of the 
Lyapunov functional H{t). Let us discuss the functional 

p{x,t) 



H{t) 



dx p{x, t) In ■ 



Po[x) 



(C.15) 



where p{x,t) is arbitrary solution of the Fokker-Planck equation ( 1C.7I) . Using the 
normalization J_ dx p{x,t) = 1 for each t > and the fact that ln{l/y) > 1 — y fro 
y > one can write the inequality 



H{t) 



dx p{x,t) 



'j^pOM^ Poix) 



Po[x) 



p{x,t) 



>0 



(C.16) 



where the equality arises for p{x,t) = Po{x). We shall show that dH/dt < for each 
t. We take the derivative of H with respect to t and by means of the Fokker-Plank 
equation flC.7l) after some calculations we obtain the relationships 

(T^ d'^p{x,t) 



dH 

"dt 



—— = dx 



-^[X{x)pix,t)]+ 2 



^2 /.oo 

— / dxp{x,t) 



dx'^ 
' Po{x) 



In 



d 



p{x,t) dx \p{x,t) 



PJx^t) 
Po{x) _ 
Po{x 



(C.17) 
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The last integral in f lC.17p is positive for p 7^ Po ^-nd it is equal to only when 

dH dH ^ ^ 

p{x,t) = po{x). Then — — < and — — = only when p{x,t) = po{x). From (10.160 

(aiL (JjL 

dH 
and from — — < it follows that 
dt - 

Urn. p{x,t) = pq{x) (C.18) 

t—¥00 

which is exactly the essence of Observation 2. 
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